!!       rnkt2.f90
!! This is the code for serial runge-kutta iterations for 
!! the spectral code. This does not make any explicit call to
!! the fourier transform libraries, thus may go unchanged from
!! one serial machine to another. 
!! --------------------------------------------------------------
subroutine rnkt(idiag)
  use general
  use mod_2dflu
  use omplib
  use force
  use mod_part_interp
  use immersed_boundary
  implicit none
  double precision,dimension(3) :: ran
  double precision :: rk2,rk2inv,delta1,tr,tc,temph,omre,omim
  double precision :: encheck,kx,ky,n_sqr
  integer :: i1,i2,k1,k2,ksqr,ireal,iimag,j,ik,id
  integer :: ierr
  integer :: idiag
  double precision :: fomegak=0.0d0,fampk=0.0d0
  double precision, dimension(3) :: kk
  double complex, dimension(3) :: wXu,curl_wXu
!! ----------------
  factor = 2.0d0*pi/length 
!! ---------------------------------------------------------
!! first  store omega in temporary  array.  
  jac_old = omega	
!!
!!----------- 
!!
!! -evaluate the Jacobean ----------------------------
  call eval_jac(1,idiag) 
!! -----------------------------------------------------------
  n_sqr = dfloat(n1*n2)
  encheck = 0.0d0
!
!DM temporarily commented out
!  if (lfourier_forcing) fampk=famp/scale
!  this should go into the fourier_forcing subroutine
!
!$OMP PARALLEL DO PRIVATE(k2,i1,k1,ireal,iimag,ksqr,rk2,kx,ky,kk,temph,tr,tc,wXu,curl_wXu,fomegak,omre,omim)  SHARED(omega)
  do i2 = 1,n2
     k2 = (i2-1) - n2*(i2/(n1hf+1))	
     do i1 =1,n1hf
        k1 = i1 -1
        ireal = 2*i1-1
        iimag=2*i1
        ksqr = k1*k1 + k2*k2 
        if((ksqr.gt.kasqr).or.(ksqr.eq.0))then
           omega(ireal,i2) = 0.0d0
           omega(iimag,i2) = 0.0d0 
        else
           rk2 = factor*factor*dfloat(ksqr)
           temph = time_increment(ksqr)
!! -------------
           kx = factor*dfloat(k1)
           ky = factor*dfloat(k2)
           kk(1) = kx; kk(2)= ky; kk(3) = 0.
           wXu(1) = dcmplx(ukx(ireal,i2),ukx(iimag,i2))
           wXu(2) = dcmplx(uky(ireal,i2),uky(iimag,i2))
           wXu(3) = dcmplx(0.,0.)
           curl_wXu = curlF(kk,wXu)
           tr = real(curl_wXu(3))
           tc = aimag(curl_wXu(3))
!           write(*,*) 'zero1',curl_wXu(1)
!           write(*,*) 'zero2',curl_wXu(2)
!           write(*,*) 'zero2',curl_wXu(3)
!--------- Fourier space force is added here ------------
           if (lfourier_forcing) then
             call add_cont_force_fourier(k1,k2,fomegak)
             tr=tr+fomegak
           endif
!---------------------------------------------------------
           omre = jac_old(ireal,i2)
           omim = jac_old(iimag,i2)
           !! ---time stepping  -----------------------------------
           !! ------------------------------------------------------------
           omega(ireal,i2)=temph*(omre+(delta/2.0d0)*tr ) 
           omega(iimag,i2)=temph*(omim+(delta/2.0d0)*tc ) 
           !! ------------------------
        endif
     enddo
  enddo
!$OMP END PARALLEL DO
!! -the second step ----------------------------------
  !! -------evaluating the non-linear part ----------------
  call eval_jac(2,idiag,delta)
!! --the second time-stepping -------------------------------------
  encheck = 0.0d0
!
!$OMP PARALLEL DO PRIVATE(k2,i1,k1,ireal,iimag,ksqr,rk2,kx,ky,temph,tr,tc,fomegak,omre,omim,kk,wXu,curl_wXu)  SHARED(omega,ukx,uky)
  do i2 = 1,n2
    k2 = (i2-1) - n2*(i2/(n1hf+1))	
    do i1 =1,n1hf
      k1 = i1 -1
      ksqr = k1*k1 + k2*k2 
      rk2 = factor*factor*dfloat(ksqr)
      ireal = 2*i1-1
      iimag=2*i1
      if((ksqr.gt.kasqr).or.(ksqr.eq.0))then
        omega(ireal,i2) = 0.0d0
        omega(iimag,i2) = 0.0d0 
      else
        rk2 = factor*factor*dfloat(ksqr)
        kx = factor*dfloat(k1)
        ky = factor*dfloat(k2)
        temph = time_increment(ksqr)
        kk(1) = kx; kk(2)= ky; kk(3) = 0.
        wXu(1) = dcmplx(ukx(ireal,i2),ukx(iimag,i2))
        wXu(2) = dcmplx(uky(ireal,i2),uky(iimag,i2))
        wXu(3) = dcmplx(0.,0.)
        curl_wXu = curlF(kk,wXu)
        tr = real(curl_wXu(3))
        tc = aimag(curl_wXu(3))
!           write(*,*) 'zero1',curl_wXu(1)
!           write(*,*) 'zero2',curl_wXu(2)
!           write(*,*) 'zero2',curl_wXu(3)
!--------- Fourier space force is added here ------------
           if (lfourier_forcing) then
             call add_cont_force_fourier(k1,k2,fomegak)
             tr=tr+fomegak
           endif
!---------------------------------------------------------
        omre = jac_old(ireal,i2)
        omim = jac_old(iimag,i2)
!! ---time stepping  -----------------------------------
        omega(ireal,i2) =  temph*(temph*omre+delta*tr) 
        omega(iimag,i2) =  temph*(temph*omim+delta*tc)
!! ---------------------------------------------------
      endif
    enddo
  enddo
!$OMP END PARALLEL DO
!! Omega is now updated. Update the immersed boundary points and the particles here. 
  if (limmersed_boundary) call evolve_immersed_boundary(delta)
  if (lparticle) call evolve_particles(delta)
!!------------------------
end subroutine rnkt
!!****************************************************

